#include "math.h"
#include "mex.h"  
#include "matrix.h"

#define himg0(i,j,k) himg0[(i+j*ix+k*ix*iy)] 
#define himg(i,j,k) himg[(i+j*ix+k*ix*iy)] 
#define img0(i,j,k) img0[(i+j*ix+k*ix*iy)] 
#define img(i,j,k) img[(i+j*ix+k*ix*iy)]
#define newimg(i,j,k) newimg[(i+j*ix+k*ix*iy)]

#define newimg0(i,j) newimg0[i+j*ix]
#define newimg1(i,j) newimg1[i+j*ix]
#define newimg2(i,j) newimg2[i+j*ix]


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
               
    double* img0 = mxGetPr(prhs[0]);
    double* img = mxGetPr(prhs[1]);
    double* himg0 = mxGetPr(prhs[2]);
    double* himg = mxGetPr(prhs[3]);

    const int *sz = mxGetDimensions(prhs[0]); 
    int ix = sz[0];
    int iy = sz[1];
    int iz = sz[2];

    double* cc = mxGetPr(prhs[4]);

    double cr1 = cc[0];
    double cr2 = cc[1];
    double cg1 = cc[2];
    double cg2 = cc[3];
    double cb1 = cc[4];
    double cb2 = cc[5];

    double ch1 = cc[6];
    double ch2 = cc[7];
    double cs1 = cc[8];
    double cs2 = cc[9];
    double cv1 = cc[10];
    double cv2 = cc[11];

    int x, y;
    double r0, g0, b0, r, g, b, h0, s0, v0, h, s, v;

    mxArray* result1 = mxCreateDoubleMatrix(ix,iy,mxREAL);
    mxArray* result2 = mxCreateDoubleMatrix(ix,iy,mxREAL);
    mxArray* result3 = mxCreateDoubleMatrix(ix,iy,mxREAL);

    double* newimg0 = mxGetPr(result1);
    double* newimg1 = mxGetPr(result2);
    double* newimg2 = mxGetPr(result3);

    int bmax = 0;

    for(x = 0; x < ix; x++)
    {   
        for(y = 0; y < iy; y++)
        {                       
                
            r = img(x,y,0);        
            g = img(x,y,1);        
            b = img(x,y,2);                
        
            r0 = img0(x,y,0);        
            g0 = img0(x,y,1);
            b0 = img0(x,y,2);                     
        
            h = himg(x,y,0);
            s = himg(x,y,1);
            v = himg(x,y,2);        
        
            h0 = himg0(x,y,0);
            s0 = himg0(x,y,1);
            v0 = himg0(x,y,2);        
       
            newimg0(x,y) = img(x,y,0);
            newimg1(x,y) = img(x,y,1);
            newimg2(x,y) = img(x,y,2);
                    
       
            if(r > (r0*cr1) && r < (r0*cr2) && g > (g0*cg1) && g < (g0*cg2) && b > (b0*cb1) && b < (b0*cb2))
            {            
                newimg0(x,y) = 0;
                newimg1(x,y) = 0;
                newimg2(x,y) = 0;
            }                               
        
            if(h > (h0*ch1) && h < (h0*ch2) && s > (s0*cs1) && s < (s0*cs2) && v > (v0*cv1) && v < (v0*cv2))
            {            
                newimg0(x,y) = 0;
                newimg1(x,y) = 0;
                newimg2(x,y) = 0;
            }              
                
        }
    }
        
    plhs[0] = result1;
    plhs[1] = result2;
    plhs[2] = result3;

            
    return;
}